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Abstract 


We consider the solution of scattering problems for the wave equation using 
approximate boundary conditions at artificial boundaries. These conditions are 
explicitly viewed as approximations to an exact boundary condition satisfied by 
the solution on the unbounded domain. We study both the short and long time 
behavior of the error. It is proved that, in two space dimensions, no local in time, 
constant coefficient boundary operator can lead to accurate results uniformly in 
time for the class of problems we consider. A variable coefficient operator is de- 
veloped which attains better accuracy (uniformly in time) than is possible with 
constant coefficient approximations. The theory is illustrated by numerical exam- 
ples. We also analyze the proposed boundary conditions using energy methods, 
leading to asymptotically correct error bounds. 


1 Introduction 

Problems posed on unbounded spatial domains arise naturally in the study of wave 
propagation. The standard approach to the numerical solution of such problems is to 
introduce an artificial boundary and apply ‘appropriate’ boundary conditions. A vast 
literature has appeared in the past 15 years, devoted primarily to the derivation of 
boundary conditions for the wave equation. (See [9] for a recent review.) Surprisingly, 
one finds very little precise error analysis in this body of work. Notable exceptions to 
this are the paper of Halpern and Rauch [17], where error estimates in terms of the 
reflection coefficient are given, and the early work of Bayliss and Turkel [3], who give 
estimates based on the size of the computational domain. In both cases the error analysis 
is (generally) not uniform in time. In the first case, this is manifested in the requirement 
that the solution not be too smooth, while in the second it is related to the long time 
breakdown of the progressive wave expansion. 

In this work we consider both the short and long time error for finite domain approx- 
imations to limiting amplitude problems for the wave equation in exterior domains. By 
limiting amplitude problems we mean cases where the forcing becomes steady or time- 
periodic as t — > oo, and one wishes to accurately compute the solution as it approaches a 
steady or time-periodic state. (In contrast, Engquist and Halpern [7], [8] consider mainly 
the rate of decay of the approximate solution to steady-state and are not interested in 
accurately reproducing the transient behavior of the actual solution.) We concentrate on 
what turns out to be the most difficult case; zero frequency (asymptotically steady) and 
two space dimensions. However, we do give some discussion of the other problems. 
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In Section 2 we develop an asymptotic error analysis for general local in time constant 
coefficient boundary conditions at a circular artificial boundary in two space dimensions. 
(We note that most of the boundary conditions which have been proposed fall into the 
category of local in time with constant coefficients.) Conditions are analyzed as approx- 
imations to the exact boundary operator, which can be conveniently expressed at such 
boundaries. It is seen that the behavior of the exact operator for low frequencies can- 
not be approximated by operators in the class considered. This immediately leads to 
lower bounds on the error of the form for long times. These are independent of the 
coefficients in the boundary condition and of the size of the computational domain. 

In Section 3 we show how to add variable coefficient corrections to the boundary 
conditions which allow us to beat the lower bounds derived earlier. We present the 
modified conditions only for the case of a circular artificial boundary, though we do make 
some remarks in the appendix on their extension to general boundaries. All rigorous error 
analyses are given for the case of a circular inner boundary. We believe the basic results 
will hold for an arbitrary star-shaped scatterer, but we do not have complete proofs. 
These results are illustrated by numerical experiments in Section 4. In particular, for a 
very simple test case, we find that for all of the standard conditions tried, the long time 
error exceeds 5%. The use of the corrected condition results in an order ol magnitude 
improvement . 

Section 5 is devoted to the analysis of the proposed conditions lor computations in 
the annular region 1 < r < R. Using a somewhat involved sequence of energy estimates, 
we prove that for long times the error is of the order fZ _ U 2 (ln f?) 2 (ln f)~ 3 . This shows 
that the error does in fact decay more rapidly than the (lnt) -1 lower bound of the 
constant coefficient case. Moreover, we show that the major source of error resides in 
the radially symmetric component. If this component is subtracted out then the error 
becomes 0((ln R) 1/2 In t • t~ 5/2 ). Short time error estimates follow from the arguments 
of [3]. (For an energy analysis of the ‘standard’ boundary conditions see Ha-Duong 
and Joly [12].) We note that the method could be extended to improve the long time 
error estimates by simply using a better approximation to the function H defined in the 
appendix. 

In our work we normalize the wave equation c~ 2 = V 2 u to have c — 1. To obtain 
the general case one needs only replace t by ct. It is interesting to note that in the 
electromagnetic problems at moderate frequencies c is very large (of the order 10 ) so 
that our ‘large V estimates essentially always hold and the improvement from (In t) to 
(lnt) -3 is very significant. 

For more complex systems with propagating solutions, for example the compressible 
Euler equations or shallow water wave equations, the construction and mathematical 
analysis of asymptotic boundary conditions is much less developed, and we are unaware 
of error analyses even in the linear case. Nonetheless, we may expect the same difficulties 
to be present for those equations which contain the wave equation as a special case. 1 he 
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Boundary Condition 

Cauchy 

/ -♦ foo{x,y) 

/ * foo{ x< iV) ‘ sin (w£ 4>) 

Dirichlet or Mixed 

0 

I 

0(i h) 


Neumann or Impedance 

_£M 


May not converge 

°W " 


Table 1: Decay of transients in 2d from Muravei [25]. 


introduction of dissipation into the governing equations is likely to improve the long 
time estimates by enhancing the decay rates, but this effect will diminish in the small 
dissipation limit. References discussing error estimates for diffusion equations include 
[ 14 , 22 ]. 

Some of the results presented here, as well as a more general discussion of the error 
analysis of approximate boundary conditions, can be found in [13]. 

2 Asymptotic Error Analysis 

Preliminary to our study of the long time behavior of the error which results from trun- 
cating the domain, we recall the basic theorems on the long-time behavior of solutions 
in exterior domains (with convex or star-shaped boundaries). Mfe consider three distinct 
cases: the Cauchy problem with compactly supported data, asymptotically constant 
forcing, and asymptotically time-periodic forcing. In the first case the solution decays 
to zero, while in the others (under appropriate assumptions) it approaches steady or 
periodic states, the so-called principle of limiting amplitude. In three space dimensions, 
the transient behavior typically decays exponentially in time [21, 24]. In two dimensions, 
in contrast, the transient behavior can be quite persistent. In an early paper, Chen [5] 
quantified this for a plane pulse incident on a circular cylinder. General theorems on the 
decay of the transients have since been given by Muravei [25], and are summarized in 
Table 2. 

The approximation of exterior problems on bounded domains leads to further ques- 
tions concerning long time behavior. In particular: 

i. Does the approximate solution approach a steady (time-periodic) state? 

ii. What is the error in the final state? 

iii. What is the transient behavior of the error? 

An affirmative answer to the first question seems a reasonable requirement to make 
on any approximation scheme, and is closely related to the notion of dissipativity for 
boundary conditions as introduced by Barry, Bielak and MacCamy [2, 4], Concentrating 
on convergence to steady state, it is possible to make the final error small or zero by 
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choosing boundary conditions which reduce to asymptotic or exact conditions for the 
Poisson equation. (See [7, 8, 16].) It is on the third question that we focus our attention 
here. In particular we will derive lower bounds of the form ^ for the error due to 
domain truncation and the imposition of any constant coefficient, local in time boundary 
condition for a two dimensional model problem. That is, no approximation of this type 
can accurately simulate the transient behavior associated with convergence to steady 
state on domains of reasonable size. In contrast, problems in three dimensions can be 
accurately solved using such ‘standard’ boundary operators, as could two dimensional 
problems obtained by imposing axisymmetry. 


2.1 Properties of the Exact Condition and Common Approx- 
imations 


Consider the wave equation forced at the surface of a cylinder: 

d 2 u d 2 u 1 du 1 d 2 u 


dt 2 


dr 2 + r dr r 2 dO 2 ' 
du , 


r > 1, 


u(r,M) = -Qj(r, 6, 0) = 0, 


u(\,0,t) = g{6,t) = goo(0) + 0 (1) , 


oo. 


( 2 . 1 . 1 ) 

( 2 . 1 . 2 ) 

(2.1.3) 


From Muravei [25], as t — » oo, u = iXoo(r,0) + where Uqo is the solution of 

Laplace’s equation exterior to the cylinder satisfying the inhomogeneous Dirichlet con- 
dition, 1 X 00 ( 1 , 0 ) = <7oo($)’ It is convenient to express u using a Fourier series in 6 and 
Laplace transforms in i. We have: 


where 


X( m)' 


= 


9n(s), 


(2.1.4) 


( 2 . 1 . 5 ) 


K n (s) 

and K n , I n are the modified Bessel functions. 

To approximate u by the solution, u, to a problem on a bounded domain, introduce 
an artificial boundary at r = R. A general form for a local in time, constant coefficient 
boundary condition there is: 

dv 

= -b n {s)v n , (2.1.6) 

where b n is a rational function of s with real coefficients. We remark that any first order 
derivative boundary conditions which do not involve purely tangential derivatives may 
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be written in the form above. If, in addition, locality in space is required, then 6„ is also 
constrained to be a real rational function of i ■ n. Note that conditions involving higher 
order r derivatives may be put in this form by substituting (s 2 ~ Iff + £) for fe- 
We view (2.1.6) as an approximation to the exact boundary condition at r = R: 


dun 

dr 


K n (Rs) ) 


Ur, 


-b' n {s)u, 


( 2 . 1 . 7 ) 


This condition, of course, is nonlocal in space and time. The error, e — u v, is easily 
expressed in transform space. (See also [13].) It satisfies: 


d 2 e n 1 de n 

dr 2 r dr 



1 < r < R, 


( 2 . 1 . 8 ) 


en(l)=0, f^(i?) + 6 n (s)e n (i?)^ = (b n (s) - b e n (s))u n (R,s). 


Solving this yields: 

e„(r,s) = A(R,s) 


K n (rs ) I n (rs) 


[ K n (s) I n (s) J 7 K n (s) ’ 


■0n(«) 


K n (Rs) 


(2.1.9) 


( 2 . 1 . 10 ) 


where A(R, s), which is the term we control by our choice of boundary condition, is given 
by: 

A(f7,s) = ( b n (s ) - b' n (s)) ■ B(R,s), (2.1.11) 


B(R,s ) 


(b n (s) - VM) 


K n (Rs) 
K n (s) 


( sJiXRs) 
[ In(Rs) 


+ b n (s)) 


In{Rs)\ 

Us) ) 


1 


( 2 . 1 . 12 ) 


In order to analyze the approximation properties of b n (s), we need to study the 
behavior of b^(s). For short to moderate times and large domains we consider the limit 
\Rs\ ~ \z\ 1. Using standard results on the asymptotic behavior of modified Bessel 

functions (e.g. [1]) we find, as z -»oo: 


K(s) = 


\zK' n {z)_ 1/1 (4n 2 — 1) (4n 2 — 1) _ 3 

RK n (z) R\ + 2 + 8 z 8z 2 + ( ’ 


(2.1.13) 


This expansion is closely related to Friedlander’s progressive wave expansion used to 
construct boundary conditions in [3, 15]. Essentially all boundary conditions proposed 
in the literature agree with at least the first term. 

We expect, on the other hand, that the long time behavior of solutions converging 
to steady state will be governed by the behavior of the transform near s 0. (This 
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corresponds with the nonuniformity in time of the progressive wave expansion.) The 
relevant expansion of 6' is: 




JKMjV-ri + 

S + 0(z ! ), 


n = 0; 
n ^ 0. 


(2.1.14) 


Matching this expansion introduces more difficulties than matching (2.1.13). If we simply 
insist that 6 n (0) = 6* (0), we see that b n must be a function of n rather than i ■ n. This 
means that the boundary operator must be nonlocal in 8. Moreover, we cannot match 
the r^- behavior for n — 0 with a rational function of a. We shall see below that this fact 
practically eliminates the possibility of accurate long time simulations using time-local 
boundary conditions with constant coefficients. 

It is of interest to recast into the present form various boundary conditions which 
have appeared in the literature, and to study them in the large and small limits. In a 
later section we will see how these properties are reflected in their performance. 
Engquist-Majda [6], Bayliss-Turkel I [3]: 


dv dv v 
~dt + fo :+ 2R = ’ 

M») = jj( z + ^). 


. r o(z->> 

" "~1 - 1 K 1 + 0 ( ln " 

Bayliss-Turkel (2nd order) [3]: 


- -♦ oo, 
1 z) 2- 0. 


2 *d d 4j - 3 

+ 7^ + 


n 

j=i 


V |r = fi= 0, 


bn(s) = 

-H 


dt dr ’ 2 r 

/ 

i + + i* + 

R 2+1 

0(z~ 3 ) 


z — * oo, 
2 — 0 . 


Engquist-Majda (3rd order) [6]: 


d 3 v d z v 1 d z v 1 d 2 v 1 d 2 v 

+ ri.o n. n/i*j + n fti<i + ft m Q/lo 5 


drdt 2 dt 3 2 R 2 dtd8 2 2 R dt 2 2 R 3 d8 2 


(2.1.15) 

(2.1.16) 


(2.1.17) 


(2.1.18) 


(2.1.19) 

(2.1.20) 


( 2 . 1 . 21 ) 


K(s) 


1 (,3+|z 2 + 4,- 

R 2 2 



( 2 . 1 . 22 ) 


7 



Engquist-Halpern [7]: 


l _ie_| 0 ( Z ') * -» OO, 

n n 1 n3 z — > 0 

l 2 Hz* Z U ' 

(2.1.23) 

dv dv 


Tt + 8? + Kv = °’ 

(2.1.24) 

b n (j) = j(* + n), 

(2.1.25) 

6— 6 e — 1 2^ 2— »00, 

° n bn -\ 0(ln -1 z) z - 0. 

(2.1.26) 


Higdon [19] has proposed simple product boundary conditions for plane boundaries. It 
is interesting to study their naive generalization to the circular case (we emphasize that 
these are not advocated in [19]): 



(2.1.27) 


b n 



M-) 


1 (1 + c x c 2 )z 2 

R (ci + C<z)z — 1 


(1+cio-ci-ca) , , q, 

c.+cj Z T 2(c 1 +cj)» + U \ Z 

—n + 0(ln _1 z) 



(2.1.28) 

(2.1.29) 


We see that condition (2.1.18) is the most accurate in the large z limit, while (2.1.24) 
is the least. Condition (2.1.24), on the other hand, is the only one with the correct z = 0 
limit (and, hence, is nonlocal in space). This corresponds to the fact it was designed 
for steady state calculations. It is surprising that (2.1.21) matches no more terms in the 
large z expansion than (2.1.15). Note, however, that it does approximate them in the 
large n limit. Condition (2.1.27) can be made to match the first two terms for large z by 
making ci = c 2 = 1. 

Conditions (2.1.21,2.1.27) are clearly not dissipative. In fact, our expression for the 
Laplace transform of (2.1.27) displays a singularity in the right half transform plane, in- 
dicating exponential growth of the solution and, hence, the error. (This will be confirmed 
in the numerical experiments.) The singularity in the transform of (2.1.21) is at s - 0. 
This allows algebraic growth in time. For example, it may be verified that t In r satis- 
fies the wave equation, a homogeneous Dirichlet condition at r = 1, as well as (2.1.21). 
Therefore, these boundary conditions are inappropriate for long time calculations. This 
conclusion echoes that of Gustafsson [11], who considered channel-like domains. 
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2.2 Long Time Asymptotics 


All the dissipative conditions above, excepting (2.1.24), will not result in the correct 
steady state limit. Moreover, no local in time constant coefficient boundary condition 
can have 6o(a) behaving correctly near 5 = 0. By studying the s — > 0 limit of (2.1.10), we 
can develop candidate asymptotic expansions of the error as t — * oo for general rational 
approximations. Although the actual behavior of the error can be worse than these 
candidate expansions predict (as with nondissipative conditions), standard Tauberian 
theorems show that it can never be better. 

The small s behavior of each term in (2.1.10) is easily computed. We begin with the 
cases n ^ 0, which are simplest: 


e„(r,s) 


( r n_ r — n)(n _ fcn(0) ) V 

(3 + MO)) -*-(*- MO)) j 


R n (] 

which suggests, as t — ► oo, 

( r n- r -n)(n- fen (0)) 


+ 0 ( 1 ), 




.*“(« + M»» -*-“(8 -Mo)) 


j « " 


Qoc,n ~f~ ^ )• 


(2.2.30) 


(2.2.31) 


The limiting value of the error is simply the error in approximating the Poisson equation 
with the boundary conditions defined by 6 n (0). It is zero if b n ( 0) = As the function 
in parentheses is bounded by 1, we have in general that the maximum of e„(r,oo) is of 
the size R~ n and decays algebraically off the artificial boundary. Practically, this means 
that only the small n harmonics need be correctly treated for s = 0. 

For n = 0 we have: 


e 0 (r,s) 


MO) 


+ 


V 1 + f>o(0) In R In 2 
This is consistent with the t — * oo expansion: 

M0) 




lnr- ^ + 0(1). 


e 0 (r,t) 


l+b o (0)lnR 


+ H(2e ^<) ) lnr • + o(< Mr 


(2.2.32) 


(2.2.33) 


where, 'H(t) is a function whose Laplace transform behaves like (sins)" 1 as s — 0, 
0((slns) -1 ) as s — + oo and has no other singularities in the closed right half plane. It is 
given by: 

(2 - 2:,4) 

Its asymptotic behavior is: 


H{t) = 


— +^— 
lnr (lnr) 2 


+ 0((lnr) 3 ) 


1 

In (re' 1 ) 


+ 0((lnr)- 3 ), 


r » oo. 


(2.2.35) 
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(The construction of 7f and its asymptotic expansion are given in Appendix A.) 

Note the very slow decay of this error with R, r and t. In particular, it may not be 
detected by varying R, for R large. Also note that the second term in the expansion is 
independent of the boundary condition. Therefore, it will be present for any dissipative 
condition, and hence provides a lower bound on the attainable long time accuracy for 
local in time, constant coefficient operators. This fact follows from: 

Theorem 1 Let e = u - v be the error resulting from domain trucation with the bound- 
ary condition (2.1.6) and let e 0 (r,t) = ^/ 0 2ir e(r,6,t)dd. Suppose /im t _ oo e 0 (r, t) exists. 
Then: 

lim sup In t ■ |e 0 (r,t) - e 0 (r,oo)| > lnr • \ 9oofi \. (2.2.36) 

f— >oo V ' 

Proof: The proof follows from the Hardy- Littlewood-Karamata Theorem (e.g. [28, 7.1 
Thm. 3]). First of all, by the Abelian theorems for Laplace transforms, 


e 0 (r,oo) 


( MQ) \ . 

\l + b o (0) In Rj lnr 


Qoo ,0 • 


(2.2.37) 


Let w(r,t) = ±(e 0 (r,oo)-e 0 (r,t)). Then w(r,s) is given by ± the second term in (2.2.32) 
as » — v 0. Suppose (2.2.36) is false. Then for some K > 0, w = w + is eventually 
positive. Moreover, 


-/ \ In r • o + K 

w(r,s)~ , s ---> 0. 

s In s 

The Hardy- Littlewood-Karamata Theorem then implies: 


(2.2.38) 


In t 

~T 



lnr • Iflfoo.ol + K, t-*oo. 


(2.2.39) 


which implies, 


In tr l 

~ j ^(r,p)dp -4 lnr ■ | goo0 | 


(2.2.40) 


This contradicts the supposition that (2.2.36) does not hold, completing the proof. O 

Again, the obstacle to proving that (2.2.33) holds in general is the possibility of worse 
behavior, for example nondecaying oscillations or even growth of the error. We will see in 
the numerical experiments that it accurately predicts the behavior of the error for most 
of the conditions listed above. 

It is worthwhile to consider the sensitivity of these results to the details of the model 
problem. The choice of a cylindrical scatterer is of no particular importance - one can 
study the solution for more general shapes using integral equation techniques. Of course, 
if there are trapped rays leading to slower decay rates for the transient solution we may 
expect even slower decay of the error. Many authors use rectangular artificial boundaries, 
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and formulate their conditions in a way that is boundary dependent. This can have some 
effect on the behavior of the conditions. For example, condition (2.1.27) then leads to 
algebraic rather than exponential growth. However, the basic conclusion concerning the 
slow decay of the error for local in time, constant coefficient conditions remains valid. 
More details can be found in Appendix B. 


3 Construction of Improved Conditions 


In this section we construct variable coefficient boundary conditions which have better 
long time behavior than can be attained by constant coefficient operators and, in addition, 
match a number of terms in the progressive wave (far field) expansion. This guarantees 
their accuracy for short to moderate times. The subsequent numerical experiments clearly 
demonstrate their efficiency. 

We assume a circular artificial boundary, r = R, and construct conditions for each 
Fourier mode. This will result in a spatially nonlocal operator, as is required if conver- 
gence to the correct steady state is to occur. We will also construct simplified approxi- 
mations which allow small errors in the final state. 

Let v = tJ n (r, £)e m *). For v n we impose a boundary condition of the form: 

KnVn - + 6n{t) ) ^ + ^ + (2 1 r + 8n{t) ) lit + = °- ( 3 - 0 - 41 ) 


To determine S n and G n we consider separately the geometrical optics and long time 
limits. The expansion (2.1.13) corresponds to the relationships: 


d 2 u n d 2 u n 1 du n 4 n 2 - 1 „_ 3 

+ -srr + + — u n = 0{R 3 u n ), 


dtdr 


dt 2 

du _ 


2 R dt 


8ft 2 


du n 1 

dl + ~9t + 2R Ur 


0(R~ 2 u n ) 


Comparing these with (3.0.41) leads immediately to the condition: 


(3.0.42) 

(3.0.43) 


/ 1 \ 1 4n 2 — 1 

6n= \ Gn ~2Rj ~8 W~' (3 '°‘ 44) 

with the formal result: 

K n u n = 0((R~ 3 + S n R- 2 )u n ). (3.0.45) 

(This will become 0(R~ 3 u n ) as S n will scale like .ft -1 .) 

To choose G n we consider the long time - near steady state behavior. For (3.0.41) we 
have: 

3v 

+ G n (t)v n m 0, (3.0.16) 
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whereas (2.1.14) implies: 


dUn 

dr 


n 

+ 5 ttn 


o(t ! ), n t^O, 


( 3 . 0 . 47 ) 


duo . 

~dr + Jo ^ ” P) u <>( R ’P) d P = °( rl )> ( 3 . 0 . 48 ) 

where Q(a) = — (iZ(ln (^) + 7)) 1 + 0(^ 2 ). Equation (3.0.47) may be approximated using 
a constant G n - 

~ n (2n + 1) 

Gn = R’ 6n= 4R ' n ^ °' ( 3 . 0 . 49 ) 

Equation (3.0.48), on the other hand, cannot be well approximated by (3.0.41) with 
constant G 0 . (This is, of course, the content of Theorem 1.) However, as u 0 (R,t) — > 
u 0 (#,oo), we have u 0 = lu 0 (i2,oo)(l + o(l)) so that: 


6 * u o = 


u 0 {R , oo) 
i?s(ln(f) + 7 ) 


(l + o(l)). 


Therefore, if we choose a decaying function G o (0 so that: 

1 


Go = — 


**(■»(?) + 7) 


(1 + o(l)), 


( 3 . 0 . 50 ) 


( 3 . 0 . 51 ) 


then CtoUo w iU agree with Q * u 0 to leading order for small s. We have already constructed 
a function with such behavor, 7i. We avoid using 7i directly, as we have no fast means 
of evaluating it, but instead use its asymptotic expansion. That is, we take: 


Go(t) ~ 




+ 0((ln<)- 3 ). 


( 3 . 0 . 52 ) 


(Note we have absorbed into one term the first two terms in the expansion of 7 i in 
inverse powers of the log of its argument. This order of agreement determines the form 
of our final error estimate. A more accurate approximation would produce a more rapid 
asymptotic error decay.) Keeping in mind that Go — (2/2) 1 ^ 0 must be enforced, we 
finally have: 


G °(‘)=Jt7. < = 2 D> 


R 


R ln( 

M0 = ^ 


R 


4 R 1 — 2(lnC)- 1 ' 

(Here D is a free parameter. In the experiments it is taken to be 4.5fl. 
Putting this all together we have our final form. 


( 3 . 0 . 53 ) 

( 3 . 0 . 54 ) 
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Uniform Asymptotic Boundary Condition 


' d , A (dv dv v s 
dl + m {Fr + al + 2R, 


1 d 2 v v 


2R 2 dB 2 8 R 2 


0 , 


(3.0.55) 


where the nonlocal (in 9) time-dependent operator V is defined by: 


\n=0 


where 


w(9,t) = 3 ? ( ^2 w n (t)e 


*(t)e ine j , 

(3.0.56) 

e in9 ) , 

(3.0.57) 


\n=0 


and S n is defined by (3.0.49) or (3.0.54). 

The main complication in using this condition, in comparison with other second order 
operators, is the application of the nonlocal operator X?. Of course this can be done using 
fast Fourier transforms, so that actual computational work is negligible. (See also Keller 
and Givoli [20] for the implementation of nonlocal conditions.) According to our analysis 
of the steady state error, the accurate approximation of 6* (0) is only import ant for small 
n. This may be used to simplify V. A first (local) approximation is to replace V by 
S 0 . The dominant part of the steady state error will then correspond to n — 1, have a 
maximum of 0(R~ 1 ) and decay like The first two terms may be incorporated in an 
operator D\ : 

Viw(0,t) = So(t)w(t) + S\(w(9,t) — w(t)), u>(t) = — f w($,t)d9. (3.0.58) 

Z7T JO 


Here we expect a maximum error of 0(i?~ 2 ) decaying like (^) 2 * Of course, this approach 
can be generalized to an arbitrary finite number of modes. We emphasize, however, there 
are no substantial savings in computational work for circular boundaries, though there 
may be in the general case. 

In Section 5 we prove the well-posedness of the resulting initial-boundary value 
problem using energy estimates, and also prove that the transient error decays faster, 
0((lnf)“ 3 ), than the lower bound established above for constant coefficient operators. A 
similar construction for three dimensional exterior problems is given in Appendix C. 


4 Numerical Experiments 

We have carried out a number of simple numerical experiments to test and illustrate the 
results derived in the preceding sections. These are necessary as our long time analysis 
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is asymptotic and doesn’t provide precise error bounds. We solved individually problems 
for various Fourier coefficients: 


In 


d 2 v n 1 d ( dv n \ n 2 

d t* rdr V* dr ) r 2 ^’ 1 < r < R ' 

(4.0.59) 

Ov 

* n M) = -^M) = °, V n (l,t) = g(t). 

(4.0.60) 

most of the experiments we have taken, for some u > : 


, cos27ru>t 

9(t) = 1 

1 + t 2 

(4.0.61) 


All simulations were carried out to t = 500, which was 50000 time steps. A standard 
second order, centered finite difference discretization was employed with a uniform spatial 
mesh and £ = .5. The solution is well-resolved by this mesh, particularly for long times 
when the solution becomes very smooth. The difference equations are: 


(At) 2 ~ 2v »( r n*) + v n (r { ,t - At)) =-- 

(Ar) 2 ( r/ Vn ( r<+1 ’*) - 2v n( r »>t) + > f )) - ^ v n (ri,t). (4.0. 


62) 


The boundary conditions used in the comparisons are (2.1.15), (2.1.18), (2.1.24) and 
(3.0.41). We also display a single simulation using (2.1.27) to illustrate the blow-up of 
the solution as f -» oo. For reference we list in Appendix D the differencing used for 
each of these, all formulas being accurate to second order. For most of the results shown 
we have taken R = 8. Note the we have carried out a large number of experiments, 
using other boundary conditions and varying R and other parameters. What we present 
is both qualitatively and quantitatively representative of what was generally observed. 

The advantage of solving the modal equations is the opportunity it provides for very 
long time and (for comparison) very large domain simulations at low computational cost. 
In particular, we have solved the difference equations with R = 252 and the same uniform 
mesh (12551 points). The error due to the boundary condition is approximated by 
comparing these solutions every 20 time steps at the point r = 1.8. Of course, the errors 
we compute for each mode may be superposed to give the errors for a full simulation. 
We are planning the development of a true multidimensional code for general star-shaped 
bodies, implementing the boundary conditions given here. 

Figures 1-3 display the errors for n = 0, which is the most difficult case. The behavior 
of the error for short times is clearly consistent with the degree to which the expansion 
(2.1.13) is matched by the boundary condition. That is, (2.1.18) is best for short times 
followed by (3.0.41), (2.1.15) and (2.1.24). These differences are somewhat accentuated 
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Figure 5: Asymptotic Expansion vs. True Error, Condition (2.1.24). 


as the frequency, u>, is made larger. The long time error for all the constant coefficient 
conditions is seen to be fairly large, ranging from 4% to 8% at t — 500. (Note that the 
error from (2.1.24) is still largest, though its limiting value is 0.) The error resulting 
from the corrected conditions, (3.0.41), is smaller than .4%, an improvement of more 
than a factor of 10. The insensitivity of these results to boundary location is illustrated 
by Figure 4, where we have doubled the domain radius. 

A comparsion of the numerically computed error and the asymptotic expansion of t he 
error (2.2.33) is presented in Figures 5-7. Plotted are the errors from Figure 1 versus: 


MQ) 

1 + fco(0) In R 



(4.0.03) 


Here, i> 0 (0) is taken from the formulas in Section 2, R = 8, r = 1.8. The expression above 
is expected to be correct to 0((lni) -3 ). We see that the candidate expansion does, in 
fact, correctly predict the behavior of the error for these particular boundary conditions. 
(Recall that in Section 2 we only prove that the error cannot decay faster.) 

Errors for n ^ 0, in contrast, are much smaller and more rapidly decaying. This is 
illustrated in Figures 8-9, corresponding to n = 1,2. Again, the overall performance of 
(3.0.41) is best. 

We have also carried out simulations for a periodic forcing, g = sin 2i xt. The results 
are shown in Figure 10. Here, the decay of the transients and the transient part of 
the error is more rapid, resulting in much smaller errors for moderate frequencies. In 
this case we find that condition (2.1.18) is most accurate, corresponding to the nearly 
invisible curve at the bottom of the figure. It is followed by (3.0.41) and (2.1.15), which 
give nearly identical results. Again, we emphasize that the absolute size of the errors is 
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Figure 10: Errors for n = 0, g = sin 2irt, R = 8. 

small compared with the earlier experiments. We would expect this to break down for 
very low frequencies, as the decay rate is 0((ut In 2 t) _1 ). 

In our final example, Figure 11, we trace the solution at a particular point computed 
with the nondissipative condition, (2.1.27). We find, as expected, that the solution 
eventually grows exponentially in time, rendering the results meaningless. 


5 Energy Estimates and Error Decay for Variable 
Coefficient Conditions 


In this section we derive energy estimates for our proposed boundary conditions. First 
of all, these imply the well-posedness of the resulting initial- boundary value problem. 
Secondly, we use them to prove that our condition is dissipative, that is that the solution 
approaches the correct steady state, and that the long time behavior of the error is 
better than can be obtained with constant coefficient operators (see Thm. 1). The 
proof is somewhat complex, relying on different techniques for the variable and constant 
coefficient parts of the boundary. To facilitate the decomposition of the error into 9- 
independent and 0-dependent pieces, we only consider the case of a circular scatterer. 

We seek to approximate the solution, u , to: 


d 2 u 

~di* 


V 2 u + /, 


(5.0.64) 


exterior to a disk, Q, of radius 1. 


We also have initial and boundary conditions (using 
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Solution 



Figure 11: Solution at r = 1.8 using (2.1.27), n = 0, u> = 1, i? = 8 . 

the usual polar coordinates r and 9 with origin at the disk center): 

du 

u — g, r = 1, u(£, 0) = w 0 (x), — (x,0) = wi(z). (5.0.65) 

We further suppose that /, Wi have compact support, that /, g approach limits as t — oo 
sufficiently fast and that all the data is smooth. 

Introducing a circular artificial boundary at r = R such that the support of the data 
is contained within the truncated domain, we approximate u by v, where v satisfies the 
wave equation with the same data as u and, in addition, (3.0.55) at r = R. The error, 
e = u — v, then satisfies: 

^ = V’e, e(£,0) = ^(i,0) = 0, xsT, (5.0.66) 

e = 0, r = 1, Be = Bu = t), r = R, (5.0.67) 

where B is the operator in (3.0.55) and T is the annulus given by 1 < r < R. We 
decompose e into two orthogonal parts, both of which satisfy the homogeneous wave 
equation with zero data away from the artificial boundary: 

e(r,9,t) = e 0 (r,t ) + e(r,0,t), (5.0.68) 


where e 0 (r, t) = (1/2 7r) / 0 2ir e(r,9,t)d9. The boundary conditions satisfied by e 0 and e are: 


v r ( deo 5^0 

al + 5 °(0 ) ( -5“ + “57 + 


>dt 


Co \ Co 


6r dt 2 R 8R 2 


= #o(0 = T r*(0,t)dO, (5.0.69) 

Z7T 70 
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de e \ 

dt + 2R) 


1 d 2 e 

2 & ae 2 


e 

8R? 


= »(«,*) = *(M)-*o(<M). 


(5.0.70) 


Here, the constant coefficient nonlocal operator V is given by the sum in (3.0.56) exclud- 
ing the n = 0 term. 

The estimates for e<j and e will be made separately, using quite different arguments. 
We note that both energy estimates could be made with a general star-shaped scatterer, 
but we would be unable to make the simple decomposition of the error into its Fourier 
modes. We expect that e 0 will represent the dominant error for long times, and so we 
will consider it first. 


5.1 Estimates for e 0 

We employ a one-parameter family of energies adapted from Muravei [26]: 


1 r R 


(P + r 2 )((^) 2 + (^Y) + 4ir de ° de ° - 9e ° 


dr 


at dr + 2feo ar “ e 


rdr , (5.1.71) 


where f = t + K and K is eventually to be chosen sufficiently large. Differentiating £ 
with respect to t, using the equation satisfied by eo, integrating by parts and applying 
the boundary condition at dSl we obtain the fundamental identity: 


£ = R 
= R 


\p + *)$£ + JS((fc)i + (*£)») + feo^J r _ R - t(%n-_ 
PSg-%- + + fg*))’ + fec,$?] r=R - f( 


(5.1.72) 


We have introduced t — t - R. 

We must now use our boundary condition at r = R to estimate the energy deriva- 
tive. As a preliminary step, we first rewrite the boundary condition in integral form by 
inverting the differential operator {m + M0) : 


de 0 

~di 


+ + ^ e °~ wJ! e ~ S ' M ' l)d ’ 1e °( R ' T)dT = (5 - L73) 


and performing an integration by parts: 
de n de 


W + W + ,( ‘ )e ° + wJo = *»('>' 

Here we have introduced: 

Z((,r) = ,(,) = ±(1 _ ^-z(t,t)), 


(5.1.74) 


(5.1.75) 
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as well as the modified data, 


$ 0 (t)= [ l e-^ s ^ )dn 9o( a ) ds - (5.1.76) 

Jo 

It is convenient to state and prove here certain facts about q which will be useful later 
on: 


Lemma 1 

q{t) > 0, q‘(t) < 0, q = 0{l/{R\n(t/ R))), t -* oo, (5.1.77) 

uniformly in R, and, for any e > 0, K = K(e)R 2 can be chosen sufficiently large such 
that: 

Proof: We recall that £o(0 is given by (3.0.54) and satisfies: 

^o(i) — ^o(0 = + ^(1/ (^/^)))’ * 00 • (5.1.79) 

This implies: 

Z(t,t) < f e- <# (*)(‘-*)ds = -^-(1 - e - 4o(t)t ) < 4 R, (5.1.80) 

Jo o 0 (£) 

which immediately yields q > 0. Differentiating q and using (5.1.80) we obtain: 

?' = -W>7t Z{t,t)= ”s7p (1 -MOZtt.0) < 0- ( 5 . 1 . 81 ) 


Now, integrating by parts, 




i m 

So (t) S>(t) 




(5.1.82) 


The first term on the right-hand side is 4/2( 1 + 0(1/ In ( tj R ))) as t — > oo and the second is 
0(R/((t/ R) In 2 ( t/ R ))). Putting this into the formula for q yields q = 0 (l/(f?ln (t/fZ))). 

Finally, we consider the terms on the left-hand side of (5.1.78). Using the expansion 
for Z and introducing t = (t + D)/R and K = 0(K/R), R — + oo, we obtain, for some 
positive constants 71 , 72 independent of R and K: 


2q < 


7i 

iZlnf 


(t + K — R)q' < - 7,(1 + 


(5.1.83) 

(5.1.84) 
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(Here we use the fact that D = 0(R).) Combining these and taking K = 0(R ) suffi- 
ciently large we see that (t+K-R)q'+2q is negative until i = 0(K / In K). Its maximum 
then, is 0(1 /(R]n K)) = 0(l/(i21n H)), and can be made arbitrarily small by choosing 
Kj(R 2 ) sufficiently large. Since the maximum of the remaining term behaves like l/K, 
the estimate holds. This completes the proof of Lemma 1 O. 

Solving for in (5.1.73,5.1.74), substituting in (5.1.72), and integrating in time from 
0 to T we obtain: 

“ 8 rL P ' {l ’’)*'))* ( 5 . 1 . 85 ) 

+Rl L 1 + 

+ Jiil ' (l .)*)) it + J*iRe c (R,t) *„(«)*. 

Our strategy from this point on is simple. We will show that the terms on the right- 
hand side of (5.1.85) involving ^ are, in aggregate, negative and we will use 8 to estimate 
the terms involving e 0 . Throughout, e and G will stand for constants independent of t, 

R » e ° 411(1 *0 sucl1 tliat e can be ma de arbitrarily small (typically at the cost of making 
K and G large). The basic estimates are presented in a sequence of lemmas: 

Lemma 2 


~2 / + < -i( f 2 q(T) + T)Re 2 Q (R,T ) (5.1.86) 

+ RhJil iRe o( R ^) dL 

Proof: Simply integrate by parts and apply the second part of Lemma 1 , choosing K 
sufficiently large. O 

Lemma 3 

L PR (jf (R,t) * (f 0 Z(<,r) ^ (fi ’ r)<ir )) dt - °- ( 5 - 1 - 87 ) 

Proof: Let, 

w ( t ) = i i ~^( R , t ), Z(t,r) = —Z(t,r). (5.1.88) 
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It is sufficient to show that: 


T t 

J w(t) J Z(t,T)w(T)d,Tdt > 0 , 


(5.1.89) 


for all functions w. By a result from the theory of Volterra equations (see Gripenberg, 
Londen and Staffans [10, 20.2 Thm. 2.2]), it is sufficient to show that, for t > r > 0, 


2>°, f <0, 


dZ_ 

dr 


> 0 , 


&z_ 

dtdr 


< 0 . 


(5.1.90) 


Recall that Z = (t + K-R)~ 2 f 0 T e~ Sl s °^ds. Clearly Z > 0 for K > R. Differentiating 
and recalling that S 0 > 0 we obtain: 


dZ 

dt 


— So{t)Z £ 0 , 


(5.1.91) 


We need 
( 5 . 1 . 80 ): 


«!£ _ , „,9z 

dtdr dr ' 

only show, then, that the r derivative of Z is nonnegative. 


(5.1.92) 


We have, using 


— = r- 3 e~£Mn)*>( r + K _ R _ 2 Z(r,r)) > r - 3 e ~ f‘ ^(^(r + K - 9 R), (5.1.93) 

which will be nonnegative for K > 9R, completing the proof. <0 
Lemma 4 


f*pR?^(R,t) •*„(()* < ^(f'el(R,T) + J* ieliR,t)dt\ (5.1.94) 

+GR 2 In R [r’iJfT’) + J* t$l(t)dt + j* p [^((]) dt j 

Proof: Integration by parts in t yields that the left-hand side of the expression above is 
equal to: 

T 2 Re 0 (R,T) • f 0 (T) - 2 [ tRe 0 (R,t) ■ %{t)dt - [ T t 2 Re 0 (R,t) ■ ^{t)dt. (5.1.95) 

Jo Jo ot 

The final result follows from the elementary product inequality ab < ea 2 -f — b 2 . O 
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Lemma 5 


[HU e-f' SoMd,, eo(R,3)dj) dt < 16.R 3 (l + J* ie 2 0 (R,t)dt. (5.1.96) 
Proof: We have, for s < t, 

(t + K)<(s + K) ^-^^- y e~r.^ < e -^. (5.1.97) 

K 

Therefore, the quantity to be estimated is bounded above by: 

(S . (< + JT)i|eo(fi, j)|) ! it, (5.1.98) 

where S(a) = (s + K)^e~vn and * denotes convolution in the time variable. By standard 
estimates for convolutions (10, 2.2 Thm. 2.2] this is bounded above by: 

^^ T S(s)d^ j* ie 2 {R,t)dt. (5.1.99) 

Integrating by parts: 

T 

[ (s + K)*e~'&ds = -4f2e~**((s + K)* + 2fi(s -f K)~^) |q 
Jo 

-4R 2 [ (s + K)-ie-™ds ( 5 . 1 . 100 ) 

Jo 

< 4RKH\ + ^). 

Substituting this into the previous expression and choosing K = 0(R\nR) sufficiently 
large yields the desired result. O 

Lemma 0 

f? tRe 0 {R,t) ■ * 0 (t)dt + ^ f? te 0 (R,t) ■ (/ 0 V /> (,,) %>( fl, s)ds) dt 

+-R 2 /o' ? ^o(0 — Si e o(A»t) + gjjy / 0 ‘ e ~ io ^ dr, e 0 (R, s)ds^j dt (5.1.101) 

< (| + sfg) J?lel(R,t)dt + GRHnRtf t*g(t)dt. 

Proof: The terms involving $0 may all be estimated using the product inequality and 
Lemma 5. Expanding the remaining terms yields: 

\ J* ie 2 0 (R , t)dt + -^J*i (jf‘ e- /•' toW'eoiR, »)da) ’ dt. (5.1.102) 

The final estimate is then a direct consequence of Lemma 5. O 
Putting all these inequalities into (5.1.85) we finally obtain: 
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Theorem 2 For any e > 0 and K = 0(R 2 ) sufficiently large there exists a constant , G, 
independent ofT, R, and ito, such that: 

e ^ + i dt + Io i2R {~d^ R,t ^j dt 

< sfgC T 2 el(R,T ) + e iel(R,t)dt) + (5.1.103) 

GR 2 In iZ (t 2 * 2 (T) + / 0 T f**(i)<« + /o T ? (^(O)* <&) 

(Here we have introduced T = T + if.) 

We now use this energy inequality to finally estimate e<> in terms of For economy 
of expression we introduce £7(T;$o) to denote the terms in parentheses on the last line 
of (5.1.103). The first step is to prove a Poincare type inequality relating e () at R to £. 


Lemma 7 For K = O(iZlniZ) sufficiently large there exists a constant L independent of 
R, T and eo such that: 

Pel(R,t) < Tin R • S(t', K). 

Proof: We begin by estimating eo in terms of f r ■ We have: 


el(R,t) = 


t R de 0( . 

J, a7 (r ’ )dr 


< In R 


i: (>•«)' 


rdr. 


By the same techniques we find that: 


rR 


J e^rdr < 


R 2 In R 


f R fde oV 

Ji \ dr ) 


rdr. 


(5.1.104) 


(5.1.105) 


(5.1.106) 


Combining this with the elementary inequalities, 


|4rf 


deo deo j 

~fr~dt' 


'deo 


de 


de 0 . 


<2rf(|^| +|^1 |, |2fe„^| < 4e 2 + ^ 


P ( de 0 


dr 


dt 


4 \ dt 


, (5.1.107) 


yields: 


£ > 


If 

2 J i 




dr 


5 

2 


dr , 


rdr. 


(5.1.108) 


Choosing K = 0(/Zln J?) sufficiently large keeps ( t — r) 2 > MR 2 In R for any M so that 
for some constant L independent of R and T we have: 


i:-a) 


rdr < L • £{t\ K). 


(5.1.109) 
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Substituting this into (5.1.105) yields the desired result. O 

Remark: The calculation used in Lemma 7 shows that, for K sufficiently large, the 

quantities jfr (%&( r ,t)) 2 dr, f*r (&(r,t)) 2 dr and |e 0 (r,t)| are all bounded by 
£(T,K) • T~ 2 . 

We are now able to estimate £ directly in terms of 'P 0 : 

Theorem 3 For any e > 0, K = K(e)R 2 can be chosen sufficiently large such that: 

£(T] K) < G(e)R\nRT'J^ i~ f ~(t;^ 0 )dt. (5.1.110) 

Proof: Making e in (5.1.103) small and using Lemma 7 to estimate e 0 we obtain: 

£(T; K)<ej* ~.£(f, K)dt -f <?(e)flln R ■ Q(T ; ¥ 0 ). (5.1.111) 

The final result follows from an application of Gronwall’s inequality. O The estimation 
of £ is completed in Theorem 8 following the estimation of if'. 


5.2 Estimates for e 


As the boundary condition for e has constant coefficients in time and the angular variable, 
we consider its Laplace transform in t and Fourier transform in 0, e n . Introducing 
differential operators, 


in 


dr 3 ^ r dr 





C n 


we have the problems (n = 1, . . . , oo): 


(4 n 2 - 1) 
8R 2 ( S + 6 n )' 


(5.2.112) 


Problem 1 

L n e n (r,s) = 0, 1 < r < R, (5.2.113) 

e n = 0, r = l; (D + C n )e n = r = R. (5.2.114) 

Our error estimates foUow in three steps. We first note that the well-posedness of the 
problem for e n follows from general principles (e.g. [27]). This implies the solution may 
be represented by its Laplace transform for fts (real part of s) sufficiently large. Second, 
as the solution of a problem with analytic coefficients (with s in the closed right half 
plane) e n must be a meromorphic function of s if 4f n is. We will first prove that there 
are no poles in the closed right half-plane, which by the preceding comment implies that 
e n is analytic. We then push the inversion contour to the imaginary s-axis, where we 
estimate the transform in terms of (Here we simply assume that the behavior at 
infinity allows this.) Finally, PlancherePs Theorem is used to give the final estimate. 
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Theorem 4 For > 0 there are no nontrivial solutions of Problem 1 with = 0. 

Proof: Suppose q is a solution of Problem 1 with zero data. Let s = ( + iq, ( > 0. 
Multiplying the equation by q and integrating from 1 to R yields: 

I n = f*r \^\ 2 dr + s 2 f*r\q\ 2 dr + n 2 f* l -\q\ 2 dr + R{s + C n (s))|q(fl)| 2 = 0. (5.2.115) 

We first note that, 


9i(s + C„(s)) — £ + 


4n 2 — 1 ( + S n 


SR 2 (t + 6 n y + r) 2 


> 0 , 


since S„ > 0. Therefore, 


0 = Rln > U 2 -V 2 + ~jp) r kl 2<ir ’ 


which implies, 






Looking now at the imaginary part, 

3/ n = 2ft j* r\q\ 2 dr + R - Z(s + Q n )\q(R ) | 2 , 

which implies, 


%(s + C n ) 


< 0. 


Writing this out and using (5.2.118) we reach a contradiction, 


S(a + Cn) _ 4n 2 — 1 

V 


I 

8 R 2 U + M 2 + V 2 - 8 - 2 


(5.2.116) 

(5.2.117) 

(5.2.118) 

(5.2.119) 

(5.2.120) 

( 5 . 2 . 121 ) 


This completes the proof o. 

We now restrict attention to the imaginary s axis where we have the identity due to 
Morawetz and Ludwig [23]: 


2 %(rDw-Lw) = V- (23? (rDtuVtu) - (s 2 |tu| 2 + |Vw| 2 ) x) (5.2.122) 
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Integrating this (with w = e n ) and using the definition of D and the boundary condition 
at r = 1 we obtain: 

l I* 

= (2» - r ((a 1 + + I^T)) If (5.2.123) 

B=(|De n |"-li^|e n p)|„ B -|f>(l)| ! . 

We now use the boundary condition to estimate the boundary term at r = R. We 
have the inequality 

2n — 1 

IG.I < -jjj- (5.2.124) 

which implies, for any /i > 0, 

|£e n | J < (1 + ^ l g "l 2 + U + (4/#*))«*. (5.2.125) 

We can use this inequality to estimate the boundary terms in (5.2.123), leading to the 
following lemma. 

Lemma 8 At r = R: 

+ ^.l ! < 5»|*„| 2 . (5.2.126) 

Proof: The lemma is a direct consequence of (5.2.125) with fi = - o. 

Substituting into (5.2.123) and applying Plancherel’s Theorem we have: 


Lemma 9 For any T > 0; 

L Ji (r* + rdr<it + ^ J 0 \ e n(R)\ 2 dt < 5R 2 n \V n \ 2 dt (5.2.127) 


Noting that e has zero average in 8, we see that Lemma 9 yields a bound on ||e||. (Here 
|| • || denotes the L 2 norm on T.) However, it can’t be used directly to get a pointwise 
bound. For that we go back to the standard energy equality (again taking account of the 
zero Dirichlet data at r = 1): 




(R,t)dt 


(5.2.128) 


The untransformed boundary condition is: 


de^ 

dr 


de n 

dt 


+ r „(*), 


(5.2.129) 
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where, 


1 4 n 2 — 1 r l 

r .W— ^O-i 


; -‘“(‘-'>e n (J?,p)dp + * n (f). (5.2.130) 


We first estimate the energy in terms of T n by substituting the boundary condition into 
(5.2.128), applying the elementary product inequality: 

1 f R ( \ @ e n ,2 , i ® e n 12 , 77,2 I 12^ J , ^ I ^ e n 12 I. ^ ^ ^ in |2 i. /» a i ni \ 


“ -L _ 

dr 1 ^ 1 3t 


+ ^W’) rdr + | /J 1^1’* < f |r.|»*. (5.2.131) 


We now estimate r„, using standard estimates for convolutions: 

f T ,n ^ J 1 /V , 2J . , (4n ! -l) ! f', 




[ T \ f e-^ t -^e n (R,p)dp\ 2 dt+ [ T \* n \ 2 dt) 
Jo Jo Jo / 


Combining these leads to an estimate which doesn’t involve the solution at the boundary. 
We introduce the fractional Sobolev norm at the boundary, IMU/2.H = (#E~ 1 ™ 3 kn| 2 ) 1/2 
and let || • || denote the L 2 norm in T. 

Theorem 5 For smooth errors e, any T > 0, and an 0(1) constant M we have the 
estimate: 


!(,r)||’ + + l;S(-.T)|* + * jf (§<*•>)'* + ~ R / 0 ' 


8 r T Ad 2 i 


< m j* 


(5.2.133) 


Proof: We begin by multiplying the inequality in Lemma 9 by ^ and adding the result 
to (5.2.131). Then, applying (5.2.132) we deduce: 

f R (.de n , de n 2 n 2 ,\ R f T ,de n . 2 4 f T f R n 4 , 

l ('ar 1 +I ^T I + 2 /o l 1t ldt+ Rj 0 l ^ TdTii 


3n 2 r T An* r 1 r'i 

+ — / \e n (R)\ 2 dt < —— j \e n (R)\ 2 dt + MRn 3 f |¥„| 2 dl. (5.2.134) 

K Jo it Jo Jo 

Subtracting the term involving e n (R) from both sides and applying Plancherel’s Theorem 

yields the desired result o. 

We can now estimate the point wise error using a Poincare inequality: 

Theorem 6 


3n 2 r T 


e 2 (r,0,T)< M\nr f ||«(t)|| 
Jo 


1 / 2 , 


(5.2.135) 


Proof: Combining (5.1.105) (with R replaced by r) with the preceding theorem imme- 
diately yields the desired result o. 
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5.3 Asymptotics of and e. 

We now use asymptotic expansions of the solution of the wave equation in combination 
with the energy estimates above to derive asymptotic error estimates. We begin with the 
well-known progressive wave expansion of Friedlander, valid for ((£ — r)/r) <C 1: 


u ~ 


Y' /j(t — r, 9) 

2-J r (2j+l)/2 ’ 


(5.3.136) 


where the functions ffip, 9) are related by: 

= J = 1,2,.... (5.3.1.17) 

All fj * s are then determined by the radiation field, /<>. Expanding in a Fourier series 
in 9 , fj = Sn /in(p)e m, , and applying the boundary condition (3.0.55) we obtain using 
(5.3.137): 

(Bu)„ ~ S„ ^ = O(r-S). (5.3.138) 

This implies, if t = R + 0(1): 

# n = 0(iH). (5.3.139) 

The short time behavior of the error can then be extracted by mimicking the scaling 
arguments of [3]. This yields the following result: 

Theorem 7 For 0 < t < (R + T), T fixed, we have 

||e|( 1( x — 0(R~ 2 ), R^oo. (5.3.140) 


(Here, || • Hij is the space-time Sobolev norm involving derivatives up to order l.) 


Estimates of the time integral of the error can then be obtained using a Poincare inequal- 
ity: 


fR+T 

/ e 2 (r,9,t)dt = 0(\nr • R~ 4 ), R — * oo, T fixed. 

Jo 


(5.3.141) 


The progressive wave expansion, however, breaks downs as t — > oo. Therefore, to 
analyze the long time behavior of the solution, we must use the long time expansion 
of Muravei [25, Sec. 5]. We decompose the solution u = u 0 + u in the same way as 
we decomposed the error, and suppose the data approaches its limit at least as fast as 
0(t ~ l ). Then we have: 


U 0 = U 0 ,oo 1 


lnr \ 
In 2 1 ) 


In r . 


(5.3.142) 
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u = Uoo{ r , 9) -f 0(ln t ■ t~ 2 ), t —> oo. 
Substituting these expressions into (3.0.55) we obtain: 

T Inf? , 


$o = f e S* 0 (tWt 

Jo 


*0 (<)_ n , Ini? 

— r> i 3 j ' ’ 


(5.3.143) 

(5.3.144) 

(5.3.145) 


(5.3.146) 


6 0 {t) ~ y R\n a t 

for the mean and, using the fact that + ^u n , «, = 0, 

$ = 0(lnt • r 3 ), 'F = 0(ln t ■ t~ 3 ). 

We also note that these are smooth functions so that the estimates can be extended to 
norms involving derivatives. Substituting these into the energy inequalities yields t he 
following error estimates: 

Theorem 8 Suppose the conditions of Theorem 3 hold and that the data approaches its 
t — > oo limit at least as fast as 0(t~ 2 ). Then for some constant L independent of R and 
T sufficiently large, 

T 2 In 3 R\ 
v R\n*T J 


S(T;K) < L 


(5.3.147) 


Proof: Substituting (5.3.145) into (5.1.103) and using the additional fact that — 

0{t~ l In -4 f ) we obtain ( L is a generic 17,T-independent constant): 

„ _ r T 2 In 2 R 
9 < L R 2 \n 6 T ’ 

(5.3.148) 

dQ ^ Tin 2 # 
dt < #Mn 6 T’ 

(5.3. 1 19) 

f T -_ t dQ ^ T 2 -‘ln 2 R 

Jo 1 dt <L R} In 6 T ‘ 

(5.3.150) 


Combining this inequality with Theorem 3 yields the desired result o. 

Lemma 7 can be used to estimate the pointwise error: 

Corollary 1 Suppose the hypotheses of Theorem 8 hold. Then, 

2/ rr ^ In 3 i? In r 
e 0 (r, i ) < L — a™ • 
ov Rln 6 T 

Finally, combining Theorem 6 with (5.3.146) we see that e decays more rapidly in 
time: 

Theorem 9 As T 


(5.3.151) 


oo, 


e 2 (r, 9, T) < In r • O 


(") 


(5.3.152) 
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A Construction and analysis of H 


We seek a function, H{t), whose Laplace transform, H(s) = (sin s)- 1 + 0(1) as s — 0, 
has no other singularities for fta > 0 and decays as fast as (sins)" 1 as s -* oo. We make 
use of the following formula, which may be found in most extensive tables of Laplace 
transforms or verified by direct computation: 


Let: 


so that: 


O f<x> r u-l+a \ 1 

p - /" — — zdu) = — - — , a > 0. 
o r(u + a) J s a lns 


f l -u~\ 

F < r > = L fn 4 ' 

JO 1 (u) 


F(s)=J-- 1 


In s s In s 

Define 7i(r ) to be the unique decaying solution of: 


dr 


-H = 


F. 


It may be represented in integral form as: 


(1.0.153) 


(1.0.154) 

(1.0.155) 


(1.0.156) 


h(t) =~r er (/o' dp - 

To compute the large r asymptotics of H, we first compute the large p asymptotics of 
F(p). By Watson’s lemma we find: 


J7(p)= b7 _ (I^)i + 0((,n '’ ) - S) - < L01,8 > 

Substituting this into the expression for Tt and integrating by parts we finally have: 

+ 0((lnr)- 3 ), r — + oo. (1.0.159) 


1 

In r 


H(t) = — — l ntn — 3>, 


(In r) J 


or/ghmal fa qe js 

OF POOR QUALITY 


B Generalizations to Noncircular Boundaries 


In this section we consider general artificial boundaries, concentrating on the construction 
of the low frequency expansion of the exact condition. Our purpose is twofold. First, we 
want to show that the lower bounds on error decay for local in time constant coefficient 
conditions cannot be improved by changing the artificial boundary shape. Second, we 
want to indicate how to construct uniformly accurate conditions for such boundaries. 
(This may be useful for scattering by bodies with large aspect ratio.) Similar ideas may 
be applied in three space dimensions. The construction of the low frequency expansion 
is discussed in more detail in [18]. 

To derive an expression for the exact condition, we study the Dirichlet problem ex- 
terior to the artificial boundary , T, which we will suppose is a level set of some smooth 
function /(x, y). After Laplace transformation we have: 


V 2 u - s 2 u = 0, 

exterior to T. By Duhamel’s principle we have: 

u(x;s) = [°° e~ st U(x,s-,t)dt 
Jo 

where U satisfies (2.0.160) and the boundary condition at T, 

U(x,s;t) = u(x,f). 

We now solve for U using integral equation techniques: 

U{x,s-,t)= f <r(y,s;t)K 0 {s\x - y\)dS y , 
v r 

leading to the singular integral equation for cr: 

f <r(y,s]t)K 0 {s\x - y\)dS y = u(x,t), x G I\ 
v r 

We recall the expansion of K 0 for small argument [1]: 


(2.0.160) 


(2.0.161) 


(2.0.162) 


(2.0.163) 


(2.0.164) 


Ko(s\x - y|) = -In \x-y\ - (3s + 0(s 2 lns), (2.0.165) 


where ft is given by: 


fi(s) = In s — In 2 + 7. 


(2.0.166) 


From [18] the small 3 asymptotics of U may be computed by substituting the small s 
asymptotics of Kq into the integral equation. We then have: 
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Theorem 10 (Hariharan and MacCamy [18]) Let fo(y) and fi(y\t) be solutions of 
the integral equations: 


J r f o(y)ln 

\x — y\dS y + Co = 0, x e r, 

(2.0.167) 


j t fo{y)dS v = 1, 

(2.0.168) 

j T S i(y;*)ln|* 

- y\dSy + Cl = -u(x,t), X G r, 

(2.0.169) 


J r fi{r>t)dSy = o. 

(2.0.170) 

Here, Co and C\ are unique constants. Let, 


W(x,s;t) = - j a 0 {y, s; t) ln |* - y\dS y - 0{s) J^a 0 (y,s\t)dS y , 

(2.0.171) 

in the exterior domain where, 



<ro(y,s\t) 

= /l at \ fo(y)- 

p(s) - Co 

(2.0.172) 

Then, 

dU dW .. 2l 2 , 

= + 0{s 2 In 2 s), s -» 0. 

On on 


U = W + 0(s 2 In 2 a), 

(2.0.173) 


To turn this into a representation of the exact boundary operator we simply compute 
the normal derivative of u, 

foo fliy 

I — dt. (2.0.174) 

Jo on 


du 

dn 


We have, 


9W = -Ma 0 = -MU ~ ST-?- M fo> 


dn ~ /3(s)-co 

where the Neumann operator M is defined by, 


Mf = *f + /(y)^- ln I* - y\ dS u 


(2.0.175) 


(2.0.176) 


Now we first note that f\ corresponds to a bounded solution to Laplace’s equation exterior 
to T with Dirichlet data — u(x,t). Therefore, 


- [°° e-'MUdt = Ku, 

Jo 


(2.0.177) 
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where K is the Dirichlet-to-Neumann map for bounded solutions of Laplace’s equation 
exterior to T. (It corresponds to the terms for the circular boundary.) The only 
dependence of the second term on u is through Cj ; 

d(i) = —Lu(x,t), (2.0.178) 


where L is some linear functional. This yields, 


Cl 

f3{s) - co 


Mf 0 = 


gg) 

M - Co 


Lu, 


(2.0.179) 


where the function Q, the constant Co, and the operator L depend only on the boundary 
and could be explicitly constructed. In transform space we therefore have for x 6 T: 

^ = Ku + Lu + 0(s 2 In 2 s), s — 1 0. (2.0.180) 

on p{&) - c 0 

Noting the dependence of /? on s, we again see a term which cannot be well approximated 
by a rational function. This means that the slow decay rates for the error which hold 
for constant coefficient conditions with circular boundaries will be present for arbitrary 
boundaries. 

Equation (2.0.180) can also be used to compute variable coefficient conditions. Again 
assuming that u approaches a steady state, 


— « Ku + Hi 2e^ c ° 7 ^)Q(x)Lu, t — > oo. 
on 


( 2 . 0 . 181 ) 


As above, the function H can be approximated by a simpler function, G(f), which has 
the same asymptotic expansion up to some number of terms. For example we could take 
G(t) = l/(co + In 2 (t + D )). 

For short times, on the other hand, we may use the geometrical optics expansion, 
which corresponds to studying the behavior of the transformed problem for s large. We 
begin with the assumption that we have chosen the boundary such that wave propagation 
is nearly normal to it. Asymptotically this would correspond to circular boundaries and 
the Friedlander expansion. But the goal of using a noncircular boundary would most 
likely be to avoid computing so far from a high aspect ratio body as required by that 
expansion. 

We further assume that /(x,y) is normalized so that |V/| “ 1 and is increasing in 
the outward normal direction. Then, by computing an expansion of the form: 

U ~ e -*/( x 'V)(^ 0 (x, y) + -g-L(x,y) +...), (2.0.182) 

s 




39 



we can compute an expansion for the normal derivative of u. Retaining as many terms 
as we did in the case of & circular boundary we obtain: 


d (du du «(®,p) \ 

m\K + Tt + ~2~ u ) 


ld 2 u 7 2 {x,y) 
2 dr 2 8 1 


0 . 


(2.0.183) 


^ denotes the tangential derivative and k, the curvature, and y 2 are given by: 


, = »v + 5 ,. = W 

dx 2 By 2 ’ \chc 2 
The uniform condition is then given by: 



«H( 


dn + dt + 2 


“) 


where T) is determined by: 


1 d 2 u 

2~d7 2 


7 a (g>y) 

8 


u — 0. 


®-H 4 £ +y )(' +0(,)fl< * ,£+ i) 


(2.0.184) 


(2.0.185) 


(2.0.186) 


(Note that the time- dependent piece represents a rank one correction, so only K + zc/2 
need be inverted.) 


C Uniform Asymptotic Conditions in 3 Dimensions 


In this section we construct conditions analagous to those constructed in Section 3 for 
the easier case of three space dimensions. Now we take a spherical artificial boundary 
and expand our approximate solution, v, in terms of spherical harmonics, Y n ( 6 , </>); 

» = «w(r f i)n. For v n we impose a condition similar to (3.0.41), but with constant 
coefficients: 


d , .> 

| dv n 

d 2 v n 

( 1 N 

l dv» 


U + \ 

f dr 

+ + 1 
dt 2 1 

\R + S 

cd| : 

i ♦ 

4* 

O 

3 

ci 

II 

P 

(3.0.187) 


This may be represented in transform space as: 


where 


-j£ + Ma)t> n = 0, 

+ 9 (5 + + S n G, 


(3.0.188) 


(3.0.189) 
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To determine the parameters 6 n , G n , we study the behavior of the exact condition, 
represented by 6* (s), for large and small s. 


bM = - 


j Q'(Rs) 


\T z 

By direct calculation we find: 


Q(Rs) ’ 


£(*) = - + -5 + 


1 n(n + 1) D _ 3 2 


i? 2i? 2 S 


+ 0(iTV 2 ), fl,s-+oo, 


Choosing: 


6 " W = TT + ° w ’ s -o. 


n + 1 n + 1 

*ri — ) ^ n — 


2fl ’ R ’ 

we have that 6 n agrees with 6* to all terms listed above. 

Using the fact that: 

+ = L ‘ Yn - ~ n ( n + 
sm 8 0<p 2 sin 0 afl 8 

we obtain a spatially nonlocal boundary condition in a fairly simple form: 

(d 


dt 


+ 


/ dv dv v\ 1 

v ) (s + Tt + r) ~ w L ‘ v = °’ 


where the spatially nonlocal operator V is defined by: 

Vw(8, <f>) = 6 n w n Y n (8, <j>), w(8, <f>) = w„V' n (^, <£). 


n=0 


n=0 


(3.0.190) 

(3.0.191) 

(3.0.192) 

(3.0.193) 

(3.0.194) 

(3.0.195) 

(3.0.196) 

(3.0.197) 


Use of this operator will result in approximations which are accurate for short to mod- 
erate times and exact at steady state. Given the exponential convergence to steady state 
for three dimensional problems, we expect the approximation to be uniformly accurate. 
However, the application of the nonlocal operator D will be significantly more expensive 
than in two dimensions, though still a small cost in comparison with the full solution pro- 
cess. (See again [20].) It is therefore useful to consider simpler approximations. Again, 
the steady state error behavior for general approximations is easy to obtain. For the 
nth harmonic the maximum error scales like and decays like ( ^)". Therefore, 
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the local operator defined by replacing 'D with do will have a steady state error which 
is 0(R~ 2 ) at the boundary and decays linearly into the interior. Replacing V with the 
simpler nonlocal operator V\ defined by: 

Viw(9,<f>) = 8 0 w + 6i(w — w), w = -j- f [ w(9, <f>) sin 9d9d<t>, (3.0.198) 

4tt J-* Jo 

yields an 0(R~ 3 ) steady state error decaying quadratically. 


D Differencing of the Boundary Conditions 


For brevity let W = v n (/Z — iAr,t + j'Af). Solution values off the grid are obtained by 
linear interpolation. 

Differencing of (2.1.15): 


X7 („*•> - »*•“) + i (»°'i - «’■*) + 


1 


At ^ 1 At 

Differencing of (2.1.18): 

1 . 1 , 

+ 77T~TT^2 + As + 


2 (R - .5Ar) 


V2' 2 = 0; 


(4.0.199) 


-(^4 + Ab) + 


: Ae = 0, (4.0.200) 


4(At) J ' 4(Ar) 2 ’ ° ' 8(J? — Ar) v (R - Ar) 2 

Ai = ((v 0>1 - 2v°'° + u 0 '" 1 ) + 2(v u - 2v 1,0 + v 1 ’- 1 ) + (v 2 ' 1 - 2v 2fi -i- n 2 " 1 )) , (4.0.201) 
A 2 = ((v 0>1 - 2V 1 ' 1 + v 2 - 1 ) + 2(v°-° - 2V 1 ' 0 + u 2 - 0 ) + (v 0 '" 1 - 2V 1 '- 1 + v 2 '" 1 )) , (4.0.202) 

(4.0.203) 




A 4 = — ((v 0,1 + 2V 1 ' 1 + v 2,1 ) — (v 0, 1 -f 2V 1, 1 4- v 2 ' a )), 
A 5 = --^-((v 0,1 + 2v 0,0 + v 0 ’- 1 ) — (v 2,1 + 2v 2,0 + v 2 ’ -1 )), 


(4.0.204) 

(4.0.205) 
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(v 0 - 1 + v 2 ' 1 + v°’~ l + V 2 ’" 1 + 2(v 1,1 + v°'° + w 1 '- 1 + v 2 -- 1 ) - yt; 1 ' 0 ) ; (4.0.206) 


Differencing of (2.1.24): 


n 


(v*’ 1 — t?i’ 0 ') H — — (v 0, » — H 

At V Ar v + (R _ ,5Ar) 


v*' j=0; 


Differencing of (3.0.41): 


Ai + A2 + A3 = 0, 


(4.0.207) 

(4.0.208) 
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* - mk; ( ( * 0 '' - vU) " + - 2 " ! '° + ” i " ) ' 


((v°-i - v 1 ' 1 ) + 2(v°’° - t, 1 - 0 ) + (v 0 '- 1 - v 1 '- 1 )) 


( 4 . 0 . 209 ) 

(4.0.210) 


A 3 = ^(«-5Ar) _ v l-i) + 6 n (t)G n (t)v?'°. (4.0.211) 

Our implementation of (2.1.27) also followed this form, noting it is equivalent (for c i = 
c 2 = 1) to (3.0.41) with G„ = 0 and 6 n = 2(fi _ 1 5Ar ) - 
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